We obtain different estimates for the standard deviations of the parameter estimates. In particular, we compare estimates under the assumption of correct model specification with robust ones.
Under correct specification, the information matrix matrix may be calculated as the
stats::optim()numDeriv::hessian()Allowing for misspecification, the robust standard errors are obtained from the covariance matrix \(\Omega^{-1}\left(\hat{\theta_T}\right) I\left(\hat{\theta_T}\right) \Omega^{-1}\left(\hat{\theta_T}\right)\). Again, the Hessian may be calculated as output of stats::optim() or in an analytic way from numDeriv::hessian().
First, we investigate the OPG and the different versions of the Hessians quantitatively and by plotting the values, see section OPG and Hessian and in particular Plot Parameters. Subsequently, we compare the diagonal elements of the OPG and the different versions of the Hessian in [Diagonal Elements]. Next, we evaluate the implied standard deviations for the system and noise parameters quantitatively and with plots in Compare Standard Deviations. Finally, we provide the results, i.e. each part of the parameter estimates together with their robust and non-robust standard errors, see Results.
Some distribution parameters are not estimated very precisely as one can see from Noise Parameters. The sub-matrix pertaining to the distribution parameters is close to singular.
Load Blanchard and Quah (1989) dataset and its dimensions.
DATASET = readRDS(params$FILEPATH_DATASET) %>%
as.matrix()
DIM_OUT = dim(DATASET)[2]
N_OBS = dim(DATASET)[1]
Load the template and the deep parameters for the best model.
PARAMS = read_rds(glue("{params$DIRPATH_MOD}bf_params_12_2.rds"))
TMPL = read_rds(glue("{params$DIRPATH_MOD}bf_tmpl_12_2.rds"))
Fill the template
params_opt_mod = fill_tmpl_whf_rev(theta = PARAMS,
tmpl = TMPL)
Print AR parameters in wide matrix format.
params_opt_mod$polm_ar %>% matrix(nrow = 2)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 -0.1073603 -0.2770627
## [2,] 0 1 0.4921382 -1.0125875
Print MA backwards parameters in wide matrix format.
params_opt_mod$polm_ma_bwd %>% matrix(nrow = 2)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0.0650218 0.74866119
## [2,] 0 1 0.2593842 0.00153858
Print MA forwards parameters in wide matrix format.
params_opt_mod$polm_ma_fwd %>% matrix(nrow = 2)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0.5610642 1.6709962
## [2,] 0 1 -0.2414728 -0.2956298
Print static matrix B
params_opt_mod$B %>% matrix(nrow = 2)
## [,1] [,2]
## [1,] 0.9357336 0.2888584
## [2,] -0.2234818 0.1819431
Print distribution parameters, \(\lambda, p, q\) for each variable (in columns)
params_opt_mod$distr %>% matrix(ncol = 2)
## [,1] [,2]
## [1,] -0.5243178 0.1492405
## [2,] 1.4924937 1.9211332
## [3,] 85616.3114036 7.8907217
Number of deep parameters parameters
TMPL$ar$n_par
## [1] 4
TMPL$ma_bwd$n_par
## [1] 4
TMPL$ma_fwd$n_par
## [1] 4
TMPL$B$n_par
## [1] 4
TMPL$distr$n_par
## [1] 6
ix_params_sys = 1:12
ix_params_sys_ar = 1:4
ix_params_sys_ma_bwd = 5:8
ix_params_sys_ma_fwd = 9:12
ix_params_noise = 13:22
ix_params_noise_b = 13:16
ix_params_noise_sgt = 17:22
Obtain likelihood function for given dataset and template.
ll_fun_sgt = ll_whf_factory(data_wide = t(DATASET),
tmpl = TMPL, shock_distr = "sgt",
use_cpp = TRUE)
Recalculate the optimisation with argument hessian = TRUE.
N_PARAMS = length(PARAMS)
out = optim(par = PARAMS,
fn = ll_fun_sgt,
method = "L-BFGS-B",
lower = c(rep(-Inf, N_PARAMS), rep(c(-1, 1.5, 2), DIM_OUT)),
upper = c(rep(Inf, N_PARAMS), rep(c(1, Inf, Inf), DIM_OUT)),
control = list(maxit = 100), hessian = TRUE)
Check output
out$par - PARAMS
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Calculate OPG
ll_grad = ll_whf_factory_grad_helper(data_wide = t(DATASET),
tmpl = TMPL)
J = jacobian(func = ll_grad,
x = PARAMS)
opg = J/sqrt(N_OBS)
opg = crossprod(opg)
if(params$ZERO_RESTR){
opg[ix_params_sys, ix_params_noise] = 0
opg[ix_params_noise, ix_params_sys] = 0
}
opg_inv = ginv(opg)
Calculate Analytic Hessian with numDeriv package.
ha = hessian(func = ll_fun_sgt,
x = PARAMS)
if(params$ZERO_RESTR){
ha[ix_params_sys, ix_params_noise] = 0
ha[ix_params_noise, ix_params_sys] = 0
}
ha_inv = ginv(ha)
Calculate Hessian via optim
hopt = out$hessian
if(params$ZERO_RESTR){
hopt[ix_params_sys, ix_params_noise] = 0
hopt[ix_params_noise, ix_params_sys] = 0
}
hopt_inv = ginv(hopt)
Hessian (inverse of information matrix under correct specification) and its normalised version.
In the first column are the
and in the second column their normalised versions.
par(mfcol = c(3,2))
opg %>% corrplot(is.corr = FALSE)
ha %>% corrplot(is.corr = FALSE)
hopt %>% corrplot(is.corr = FALSE)
opg %>% cov2cor() %>% corrplot() %>% suppressWarnings()
ha %>% cov2cor() %>% corrplot() %>% suppressWarnings()
hopt %>% cov2cor() %>% corrplot() %>% suppressWarnings()
par(mfcol = c(3,2))
opg[ix_params_sys, ix_params_sys] %>% corrplot(is.corr = FALSE)
ha[ix_params_sys, ix_params_sys] %>% corrplot(is.corr = FALSE)
hopt[ix_params_sys, ix_params_sys] %>% corrplot(is.corr = FALSE)
opg[ix_params_sys, ix_params_sys] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
ha[ix_params_sys, ix_params_sys] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
hopt[ix_params_sys, ix_params_sys] %>%cov2cor() %>% corrplot() %>% suppressWarnings()
opg[ix_params_sys, ix_params_noise] = 0
opg[ix_params_noise, ix_params_sys] = 0
opg[ix_params_noise_b, ix_params_noise_sgt] = 0
opg[ix_params_noise_sgt, ix_params_noise_b] = 0
opg[ix_params_noise_sgt[1:3], ix_params_noise_sgt[4:6]] = 0
opg[ix_params_noise_sgt[4:6], ix_params_noise_sgt[1:3]] = 0
opg %>% zapsmall()
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 8.30260 5.51231 -0.37724 2.80218 -8.21639 -4.26124 2.53252 3.27255
## [2,] 5.51231 29.31015 -3.27003 -1.30250 -6.63013 -27.34849 1.12524 7.23095
## [3,] -0.37724 -3.27003 4.37678 -8.81453 1.63817 -0.32804 -0.76047 -0.82736
## [4,] 2.80218 -1.30250 -8.81453 52.56240 -5.27292 19.03739 2.65373 -2.56562
## [5,] -8.21639 -6.63013 1.63817 -5.27292 8.75659 4.66961 -2.52694 -3.42208
## [6,] -4.26124 -27.34849 -0.32804 19.03739 4.66961 32.68539 0.36090 -7.17842
## [7,] 2.53252 1.12524 -0.76047 2.65373 -2.52694 0.36090 1.27433 1.17358
## [8,] 3.27255 7.23095 -0.82736 -2.56562 -3.42208 -7.17842 1.17358 3.77695
## [9,] 0.02605 0.19201 -0.54652 0.30418 -0.15658 -0.05488 0.18029 0.27542
## [10,] -1.47012 1.80855 -0.81566 0.50156 1.07335 -1.91806 -0.52383 0.11658
## [11,] 0.20368 0.02012 0.39924 -0.99184 -0.07523 -0.42292 -0.07180 0.01449
## [12,] -0.48979 -1.34146 0.94565 -3.17236 0.80975 -0.11514 -0.42854 -0.45038
## [13,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [14,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [15,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [16,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [17,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [18,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [19,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [20,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [21,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [22,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
## [1,] 0.02605 -1.47012 0.20368 -0.48979 0.00000 0.00000 0.00000 0.00000
## [2,] 0.19201 1.80855 0.02012 -1.34146 0.00000 0.00000 0.00000 0.00000
## [3,] -0.54652 -0.81566 0.39924 0.94565 0.00000 0.00000 0.00000 0.00000
## [4,] 0.30418 0.50156 -0.99184 -3.17236 0.00000 0.00000 0.00000 0.00000
## [5,] -0.15658 1.07335 -0.07523 0.80975 0.00000 0.00000 0.00000 0.00000
## [6,] -0.05488 -1.91806 -0.42292 -0.11514 0.00000 0.00000 0.00000 0.00000
## [7,] 0.18029 -0.52383 -0.07180 -0.42854 0.00000 0.00000 0.00000 0.00000
## [8,] 0.27542 0.11658 0.01449 -0.45038 0.00000 0.00000 0.00000 0.00000
## [9,] 2.44750 4.33006 -0.44009 -0.53731 0.00000 0.00000 0.00000 0.00000
## [10,] 4.33006 19.07369 -0.67460 -1.99645 0.00000 0.00000 0.00000 0.00000
## [11,] -0.44009 -0.67460 0.30339 0.38985 0.00000 0.00000 0.00000 0.00000
## [12,] -0.53731 -1.99645 0.38985 2.05647 0.00000 0.00000 0.00000 0.00000
## [13,] 0.00000 0.00000 0.00000 0.00000 2.76124 3.10288 1.01643 -0.62951
## [14,] 0.00000 0.00000 0.00000 0.00000 3.10288 26.07433 3.52805 -7.89943
## [15,] 0.00000 0.00000 0.00000 0.00000 1.01643 3.52805 2.29553 3.91247
## [16,] 0.00000 0.00000 0.00000 0.00000 -0.62951 -7.89943 3.91247 31.59374
## [17,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [18,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [19,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [20,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [21,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [22,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [,17] [,18] [,19] [,20] [,21] [,22]
## [1,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [2,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [3,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [4,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [5,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [6,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [7,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [8,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [9,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [10,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [11,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [12,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [13,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [14,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [15,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [16,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [17,] 0.87789 0.07510 0 0.00000 0.00000 0.00000
## [18,] 0.07510 0.13163 0 0.00000 0.00000 0.00000
## [19,] 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [20,] 0.00000 0.00000 0 0.50289 -0.00350 -0.00029
## [21,] 0.00000 0.00000 0 -0.00350 0.05719 0.00166
## [22,] 0.00000 0.00000 0 -0.00029 0.00166 0.00005
opg %>% ginv() %>% zapsmall()
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 7.74 -0.78 -2.40 0.24 6.79 -0.70 -3.40 0.28 -0.20 -0.02 -1.11 -0.42
## [2,] -0.78 2.43 0.46 -0.82 -1.73 2.68 -2.64 -0.05 0.09 0.04 -0.46 0.34
## [3,] -2.40 0.46 1.22 -0.06 -2.16 0.42 0.83 0.04 0.10 0.01 -0.01 0.17
## [4,] 0.24 -0.82 -0.06 0.33 0.61 -0.93 0.91 0.07 0.00 -0.02 0.18 -0.08
## [5,] 6.79 -1.73 -2.16 0.61 6.84 -1.86 -0.88 0.25 -0.16 -0.05 -0.41 -0.51
## [6,] -0.70 2.68 0.42 -0.93 -1.86 3.07 -3.39 0.17 0.07 0.06 -0.55 0.36
## [7,] -3.40 -2.64 0.83 0.91 -0.88 -3.39 8.91 -1.22 -0.16 0.04 1.49 -0.04
## [8,] 0.28 -0.05 0.04 0.07 0.25 0.17 -1.22 1.12 -0.02 -0.01 -0.16 0.04
## [9,] -0.20 0.09 0.10 0.00 -0.16 0.07 -0.16 -0.02 0.97 -0.20 1.21 -0.18
## [10,] -0.02 0.04 0.01 -0.02 -0.05 0.06 0.04 -0.01 -0.20 0.11 -0.18 0.10
## [11,] -1.11 -0.46 -0.01 0.18 -0.41 -0.55 1.49 -0.16 1.21 -0.18 6.85 -1.02
## [12,] -0.42 0.34 0.17 -0.08 -0.51 0.36 -0.04 0.04 -0.18 0.10 -1.02 0.87
## [13,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [14,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [15,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [16,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [17,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [18,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [19,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [22,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## [1,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [2,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [3,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [4,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [5,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [6,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [7,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [8,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [9,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [10,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [11,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [12,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [13,] 0.48 -0.01 -0.27 0.04 0.00 0.00 0 0.00 0.00 0.00
## [14,] -0.01 0.08 -0.19 0.04 0.00 0.00 0 0.00 0.00 0.00
## [15,] -0.27 -0.19 1.18 -0.20 0.00 0.00 0 0.00 0.00 0.00
## [16,] 0.04 0.04 -0.20 0.07 0.00 0.00 0 0.00 0.00 0.00
## [17,] 0.00 0.00 0.00 0.00 1.20 -0.68 0 0.00 0.00 0.00
## [18,] 0.00 0.00 0.00 0.00 -0.68 7.99 0 0.00 0.00 0.00
## [19,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [20,] 0.00 0.00 0.00 0.00 0.00 0.00 0 2.01 -1.72 63.53
## [21,] 0.00 0.00 0.00 0.00 0.00 0.00 0 -1.72 158.15 -4856.51
## [22,] 0.00 0.00 0.00 0.00 0.00 0.00 0 63.53 -4856.51 167688.11
hopt[ix_params_sys, ix_params_noise] = 0
hopt[ix_params_noise, ix_params_sys] = 0
hopt[ix_params_noise_b, ix_params_noise_sgt] = 0
hopt[ix_params_noise_sgt, ix_params_noise_b] = 0
hopt[ix_params_noise_sgt[1:3], ix_params_noise_sgt[4:6]] = 0
hopt[ix_params_noise_sgt[4:6], ix_params_noise_sgt[1:3]] = 0
hopt %>% zapsmall()
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 6.80671 1.15824 1.58211 -2.24422 -6.18947 -2.71321 1.22165
## [2,] 1.15824 24.23756 -4.78356 7.13560 -3.35882 -19.07695 0.95329
## [3,] 1.58211 -4.78356 6.88489 -14.63398 0.66565 -0.92152 -0.36172
## [4,] -2.24422 7.13560 -14.63398 68.70189 -2.18633 15.49050 0.12931
## [5,] -6.18947 -3.35882 0.66565 -2.18633 6.86046 3.24671 -1.21063
## [6,] -2.71321 -19.07695 -0.92152 15.49050 3.24671 23.47750 0.20038
## [7,] 1.22165 0.95329 -0.36172 0.12931 -1.21063 0.20038 1.20988
## [8,] 1.57283 3.63632 -0.77461 -2.12066 -1.63627 -3.87901 0.60253
## [9,] -1.80201 -1.05038 -1.76138 3.38608 1.54798 2.45213 -0.21297
## [10,] 0.04229 0.60600 -0.51997 -1.46983 1.00414 0.14682 1.22563
## [11,] 0.51867 -0.39122 0.44157 -1.27603 -0.43692 -0.06173 0.10368
## [12,] -0.31442 0.44981 0.12153 -1.27525 0.00428 -1.01057 -0.29194
## [13,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [14,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [15,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [16,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [17,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [18,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [19,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [20,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [21,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [22,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
## [1,] 1.57283 -1.80201 0.04229 0.51867 -0.31442 0.00000 0.00000 0.00000
## [2,] 3.63632 -1.05038 0.60600 -0.39122 0.44981 0.00000 0.00000 0.00000
## [3,] -0.77461 -1.76138 -0.51997 0.44157 0.12153 0.00000 0.00000 0.00000
## [4,] -2.12066 3.38608 -1.46983 -1.27603 -1.27525 0.00000 0.00000 0.00000
## [5,] -1.63627 1.54798 1.00414 -0.43692 0.00428 0.00000 0.00000 0.00000
## [6,] -3.87901 2.45213 0.14682 -0.06173 -1.01057 0.00000 0.00000 0.00000
## [7,] 0.60253 -0.21297 1.22563 0.10368 -0.29194 0.00000 0.00000 0.00000
## [8,] 2.54610 0.53657 0.80833 0.05397 0.61825 0.00000 0.00000 0.00000
## [9,] 0.53657 4.84800 4.80466 -0.90319 -0.64090 0.00000 0.00000 0.00000
## [10,] 0.80833 4.80466 21.27074 -0.58546 -3.19909 0.00000 0.00000 0.00000
## [11,] 0.05397 -0.90319 -0.58546 0.36647 0.29328 0.00000 0.00000 0.00000
## [12,] 0.61825 -0.64090 -3.19909 0.29328 1.87076 0.00000 0.00000 0.00000
## [13,] 0.00000 0.00000 0.00000 0.00000 0.00000 2.65480 0.59093 0.66053
## [14,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.59093 19.30288 3.20398
## [15,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.66053 3.20398 3.75069
## [16,] 0.00000 0.00000 0.00000 0.00000 0.00000 -1.06378 -5.10908 2.86177
## [17,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [18,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [19,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [20,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [21,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [22,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## [1,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [2,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [3,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [4,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [5,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [6,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [7,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [8,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [9,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [10,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [11,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [12,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [13,] -1.06378 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [14,] -5.10908 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [15,] 2.86177 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [16,] 31.95176 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [17,] 0.00000 2.46573 -0.14302 0 0.00000 0.00000 0.00000
## [18,] 0.00000 -0.14302 0.17574 0 0.00000 0.00000 0.00000
## [19,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [20,] 0.00000 0.00000 0.00000 0 0.44235 -0.05234 -0.00107
## [21,] 0.00000 0.00000 0.00000 0 -0.05234 0.07686 0.00238
## [22,] 0.00000 0.00000 0.00000 0 -0.00107 0.00238 0.00008
hopt %>% ginv() %>% zapsmall()
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 4.03 0.83 -1.47 -0.36 3.18 0.64 -0.93 -1.50 0.44 -0.08 0.25
## [2,] 0.83 7.33 0.06 -2.65 -2.12 8.70 -11.30 1.80 -2.22 0.58 -5.32
## [3,] -1.47 0.06 0.91 0.09 -1.25 0.15 -0.08 0.74 -0.14 0.03 -0.19
## [4,] -0.36 -2.65 0.09 1.00 0.75 -3.18 4.13 -0.61 0.84 -0.22 2.09
## [5,] 3.18 -2.12 -1.25 0.75 3.84 -2.86 3.76 -1.86 1.25 -0.33 2.67
## [6,] 0.64 8.70 0.15 -3.18 -2.86 10.49 -13.66 2.50 -2.91 0.74 -6.93
## [7,] -0.93 -11.30 -0.08 4.13 3.76 -13.66 19.06 -3.53 3.94 -1.04 8.86
## [8,] -1.50 1.80 0.74 -0.61 -1.86 2.50 -3.53 2.37 -1.33 0.22 -2.37
## [9,] 0.44 -2.22 -0.14 0.84 1.25 -2.91 3.94 -1.33 1.75 -0.39 3.89
## [10,] -0.08 0.58 0.03 -0.22 -0.33 0.74 -1.04 0.22 -0.39 0.16 -0.89
## [11,] 0.25 -5.32 -0.19 2.09 2.67 -6.93 8.86 -2.37 3.89 -0.89 13.07
## [12,] 1.00 -0.02 -0.40 -0.03 0.73 -0.08 0.33 -0.84 0.00 0.15 -1.05
## [13,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [14,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [15,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [16,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [17,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [18,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [19,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [22,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
## [1,] 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [2,] -0.02 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [3,] -0.40 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [4,] -0.03 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [5,] 0.73 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [6,] -0.08 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [7,] 0.33 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [8,] -0.84 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [9,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [10,] 0.15 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [11,] -1.05 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [12,] 1.42 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [13,] 0.00 0.41 0.01 -0.10 0.02 0.00 0.00 0 0.00 0.00
## [14,] 0.00 0.01 0.07 -0.07 0.02 0.00 0.00 0 0.00 0.00
## [15,] 0.00 -0.10 -0.07 0.39 -0.05 0.00 0.00 0 0.00 0.00
## [16,] 0.00 0.02 0.02 -0.05 0.04 0.00 0.00 0 0.00 0.00
## [17,] 0.00 0.00 0.00 0.00 0.00 0.43 0.35 0 0.00 0.00
## [18,] 0.00 0.00 0.00 0.00 0.00 0.35 5.97 0 0.00 0.00
## [19,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00
## [20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0 2.86 11.08
## [21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0 11.08 236.46
## [22,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0 -294.85 -6967.65
## [,22]
## [1,] 0.00
## [2,] 0.00
## [3,] 0.00
## [4,] 0.00
## [5,] 0.00
## [6,] 0.00
## [7,] 0.00
## [8,] 0.00
## [9,] 0.00
## [10,] 0.00
## [11,] 0.00
## [12,] 0.00
## [13,] 0.00
## [14,] 0.00
## [15,] 0.00
## [16,] 0.00
## [17,] 0.00
## [18,] 0.00
## [19,] 0.00
## [20,] -294.85
## [21,] -6967.65
## [22,] 218382.47
ha[ix_params_sys, ix_params_noise] = 0
ha[ix_params_noise, ix_params_sys] = 0
ha[ix_params_noise_b, ix_params_noise_sgt] = 0
ha[ix_params_noise_sgt, ix_params_noise_b] = 0
ha[ix_params_noise_sgt[1:3], ix_params_noise_sgt[4:6]] = 0
ha[ix_params_noise_sgt[4:6], ix_params_noise_sgt[1:3]] = 0
ha %>% zapsmall()
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 6.64779 1.38285 1.55838 0.21292 -6.13935 -2.15539 1.37702
## [2,] 1.38285 23.29748 -4.65762 4.68690 -2.95251 -18.50730 0.22495
## [3,] 1.55838 -4.65762 7.16766 -11.52276 0.54254 -0.78910 -0.36629
## [4,] 0.21292 4.68690 -11.52276 62.16895 -3.75032 14.43507 0.44700
## [5,] -6.13935 -2.95251 0.54254 -3.75032 6.84693 2.99103 -1.12235
## [6,] -2.15539 -18.50730 -0.78910 14.43507 2.99103 23.52999 0.63428
## [7,] 1.37702 0.22495 -0.36629 0.44700 -1.12235 0.63428 1.04505
## [8,] 1.57400 3.65093 -0.85760 -3.95510 -1.63601 -4.14276 0.63166
## [9,] -1.29375 -0.69725 -1.13021 -0.37457 1.34147 1.16890 -0.09826
## [10,] 0.14462 0.95949 -0.07914 -1.73853 1.06119 -0.05290 1.15974
## [11,] 0.34404 -0.55196 0.23172 -0.65132 -0.34884 0.14036 0.05165
## [12,] -0.31968 0.50839 0.02581 -1.50673 0.00026 -0.95949 -0.27621
## [13,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [14,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [15,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [16,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [17,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [18,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [19,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [20,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [21,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [22,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
## [1,] 1.57400 -1.29375 0.14462 0.34404 -0.31968 0.00000 0.00000 0.00000
## [2,] 3.65093 -0.69725 0.95949 -0.55196 0.50839 0.00000 0.00000 0.00000
## [3,] -0.85760 -1.13021 -0.07914 0.23172 0.02581 0.00000 0.00000 0.00000
## [4,] -3.95510 -0.37457 -1.73853 -0.65132 -1.50673 0.00000 0.00000 0.00000
## [5,] -1.63601 1.34147 1.06119 -0.34884 0.00026 0.00000 0.00000 0.00000
## [6,] -4.14276 1.16890 -0.05290 0.14036 -0.95949 0.00000 0.00000 0.00000
## [7,] 0.63166 -0.09826 1.15974 0.05165 -0.27621 0.00000 0.00000 0.00000
## [8,] 2.53165 0.54167 0.71112 0.05541 0.61261 0.00000 0.00000 0.00000
## [9,] 0.54167 4.18309 4.79907 -0.75394 -0.64415 0.00000 0.00000 0.00000
## [10,] 0.71112 4.79907 21.11990 -0.57324 -3.14332 0.00000 0.00000 0.00000
## [11,] 0.05541 -0.75394 -0.57324 0.34160 0.29750 0.00000 0.00000 0.00000
## [12,] 0.61261 -0.64415 -3.14332 0.29750 1.87748 0.00000 0.00000 0.00000
## [13,] 0.00000 0.00000 0.00000 0.00000 0.00000 2.61700 0.71211 0.61947
## [14,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.71211 19.33209 2.95817
## [15,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.61947 2.95817 3.69067
## [16,] 0.00000 0.00000 0.00000 0.00000 0.00000 -1.24011 -5.22227 2.89039
## [17,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [18,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [19,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [20,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [21,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [22,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## [1,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [2,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [3,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [4,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [5,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [6,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [7,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [8,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [9,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [10,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [11,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [12,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [13,] -1.24011 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [14,] -5.22227 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [15,] 2.89039 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [16,] 32.05416 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [17,] 0.00000 2.20306 -0.10606 0 0.00000 0.00000 0.00000
## [18,] 0.00000 -0.10606 0.18180 0 0.00000 0.00000 0.00000
## [19,] 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## [20,] 0.00000 0.00000 0.00000 0 0.44234 -0.05233 -0.00107
## [21,] 0.00000 0.00000 0.00000 0 -0.05233 0.07686 0.00238
## [22,] 0.00000 0.00000 0.00000 0 -0.00107 0.00238 0.00008
ha %>% ginv() %>% zapsmall()
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 6.90 -1.24 -3.02 0.13 5.92 -1.57 -0.83 -2.27 -0.05 0.07 -0.97 1.72
## [2,] -1.24 0.72 0.69 -0.15 -1.22 0.75 -0.23 0.29 0.18 -0.05 0.77 -0.43
## [3,] -3.02 0.69 1.67 -0.01 -2.54 0.81 0.31 1.21 0.20 -0.09 0.90 -0.89
## [4,] 0.13 -0.15 -0.01 0.08 0.21 -0.18 0.12 0.08 0.01 -0.01 0.05 0.02
## [5,] 5.92 -1.22 -2.54 0.21 5.44 -1.56 -0.35 -1.70 -0.05 0.02 -0.56 1.36
## [6,] -1.57 0.75 0.81 -0.18 -1.56 0.94 -0.53 0.69 -0.06 -0.01 0.21 -0.51
## [7,] -0.83 -0.23 0.31 0.12 -0.35 -0.53 2.86 -0.82 0.64 -0.17 1.00 0.20
## [8,] -2.27 0.29 1.21 0.08 -1.70 0.69 -0.82 2.57 -0.62 0.02 -0.59 -1.11
## [9,] -0.05 0.18 0.20 0.01 -0.05 -0.06 0.64 -0.62 1.11 -0.23 2.43 -0.17
## [10,] 0.07 -0.05 -0.09 -0.01 0.02 -0.01 -0.17 0.02 -0.23 0.12 -0.54 0.20
## [11,] -0.97 0.77 0.90 0.05 -0.56 0.21 1.00 -0.59 2.43 -0.54 9.68 -1.50
## [12,] 1.72 -0.43 -0.89 0.02 1.36 -0.51 0.20 -1.11 -0.17 0.20 -1.50 1.62
## [13,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [14,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [15,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [16,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [17,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [18,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [19,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [22,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## [1,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [2,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [3,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [4,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [5,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [6,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [7,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [8,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [9,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [10,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [11,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [12,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [13,] 0.41 0.01 -0.09 0.03 0.00 0.00 0 0.00 0.00 0.00
## [14,] 0.01 0.07 -0.07 0.02 0.00 0.00 0 0.00 0.00 0.00
## [15,] -0.09 -0.07 0.38 -0.05 0.00 0.00 0 0.00 0.00 0.00
## [16,] 0.03 0.02 -0.05 0.04 0.00 0.00 0 0.00 0.00 0.00
## [17,] 0.00 0.00 0.00 0.00 0.47 0.27 0 0.00 0.00 0.00
## [18,] 0.00 0.00 0.00 0.00 0.27 5.66 0 0.00 0.00 0.00
## [19,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [20,] 0.00 0.00 0.00 0.00 0.00 0.00 0 2.86 11.08 -294.85
## [21,] 0.00 0.00 0.00 0.00 0.00 0.00 0 11.08 236.46 -6967.70
## [22,] 0.00 0.00 0.00 0.00 0.00 0.00 0 -294.85 -6967.70 218384.10
c(opg[ix_params_noise, ix_params_noise],
ha[ix_params_noise, ix_params_noise],
hopt[ix_params_noise, ix_params_noise]) %>%
array(dim = c(10, 10, 3)) %>%
zapsmall()
## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 2.76124 3.10288 1.01643 -0.62951 0.00000 0.00000 0 0.00000 0.00000
## [2,] 3.10288 26.07433 3.52805 -7.89943 0.00000 0.00000 0 0.00000 0.00000
## [3,] 1.01643 3.52805 2.29553 3.91247 0.00000 0.00000 0 0.00000 0.00000
## [4,] -0.62951 -7.89943 3.91247 31.59374 0.00000 0.00000 0 0.00000 0.00000
## [5,] 0.00000 0.00000 0.00000 0.00000 0.87789 0.07510 0 0.00000 0.00000
## [6,] 0.00000 0.00000 0.00000 0.00000 0.07510 0.13163 0 0.00000 0.00000
## [7,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 0.00000 0.00000
## [8,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 0.50289 -0.00350
## [9,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 -0.00350 0.05719
## [10,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 -0.00029 0.00166
## [,10]
## [1,] 0.00000
## [2,] 0.00000
## [3,] 0.00000
## [4,] 0.00000
## [5,] 0.00000
## [6,] 0.00000
## [7,] 0.00000
## [8,] -0.00029
## [9,] 0.00166
## [10,] 0.00005
##
## , , 2
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 2.61700 0.71211 0.61947 -1.24011 0.00000 0.00000 0 0.00000
## [2,] 0.71211 19.33209 2.95817 -5.22227 0.00000 0.00000 0 0.00000
## [3,] 0.61947 2.95817 3.69067 2.89039 0.00000 0.00000 0 0.00000
## [4,] -1.24011 -5.22227 2.89039 32.05416 0.00000 0.00000 0 0.00000
## [5,] 0.00000 0.00000 0.00000 0.00000 2.20306 -0.10606 0 0.00000
## [6,] 0.00000 0.00000 0.00000 0.00000 -0.10606 0.18180 0 0.00000
## [7,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 0.00000
## [8,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 0.44234
## [9,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 -0.05233
## [10,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 -0.00107
## [,9] [,10]
## [1,] 0.00000 0.00000
## [2,] 0.00000 0.00000
## [3,] 0.00000 0.00000
## [4,] 0.00000 0.00000
## [5,] 0.00000 0.00000
## [6,] 0.00000 0.00000
## [7,] 0.00000 0.00000
## [8,] -0.05233 -0.00107
## [9,] 0.07686 0.00238
## [10,] 0.00238 0.00008
##
## , , 3
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 2.65480 0.59093 0.66053 -1.06378 0.00000 0.00000 0 0.00000
## [2,] 0.59093 19.30288 3.20398 -5.10908 0.00000 0.00000 0 0.00000
## [3,] 0.66053 3.20398 3.75069 2.86177 0.00000 0.00000 0 0.00000
## [4,] -1.06378 -5.10908 2.86177 31.95176 0.00000 0.00000 0 0.00000
## [5,] 0.00000 0.00000 0.00000 0.00000 2.46573 -0.14302 0 0.00000
## [6,] 0.00000 0.00000 0.00000 0.00000 -0.14302 0.17574 0 0.00000
## [7,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 0.00000
## [8,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 0.44235
## [9,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 -0.05234
## [10,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0 -0.00107
## [,9] [,10]
## [1,] 0.00000 0.00000
## [2,] 0.00000 0.00000
## [3,] 0.00000 0.00000
## [4,] 0.00000 0.00000
## [5,] 0.00000 0.00000
## [6,] 0.00000 0.00000
## [7,] 0.00000 0.00000
## [8,] -0.05234 -0.00107
## [9,] 0.07686 0.00238
## [10,] 0.00238 0.00008
hopt[ix_params_noise, ix_params_noise] %>% ginv() %>% zapsmall()
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0.41 0.01 -0.10 0.02 0.00 0.00 0 0.00 0.00 0.00
## [2,] 0.01 0.07 -0.07 0.02 0.00 0.00 0 0.00 0.00 0.00
## [3,] -0.10 -0.07 0.39 -0.05 0.00 0.00 0 0.00 0.00 0.00
## [4,] 0.02 0.02 -0.05 0.04 0.00 0.00 0 0.00 0.00 0.00
## [5,] 0.00 0.00 0.00 0.00 0.43 0.35 0 0.00 0.00 0.00
## [6,] 0.00 0.00 0.00 0.00 0.35 5.97 0 0.00 0.00 0.00
## [7,] 0.00 0.00 0.00 0.00 0.00 0.00 0 0.00 0.00 0.00
## [8,] 0.00 0.00 0.00 0.00 0.00 0.00 0 2.86 11.08 -294.85
## [9,] 0.00 0.00 0.00 0.00 0.00 0.00 0 11.08 236.46 -6967.65
## [10,] 0.00 0.00 0.00 0.00 0.00 0.00 0 -294.85 -6967.65 218382.47
hopt[ix_params_noise_sgt, ix_params_noise_sgt] %>% ginv() %>% zapsmall()
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.4 0.3 0.4 0.0 0.0 0.0
## [2,] 0.3 6.0 3.4 0.0 0.0 0.0
## [3,] 0.4 3.4 -2180917.9 0.0 0.0 0.0
## [4,] 0.0 0.0 0.0 2.9 11.1 -294.9
## [5,] 0.0 0.0 0.0 11.1 236.5 -6967.7
## [6,] 0.0 0.0 0.0 -294.9 -6967.7 218382.5
eigen(hopt[ix_params_noise, ix_params_noise])$values
## [1] 3.389182e+01 1.856179e+01 3.248406e+00 2.474626e+00 1.958111e+00
## [6] 4.496991e-01 1.668400e-01 6.958517e-02 4.574457e-06 -4.585225e-07
eigen(hopt[ix_params_noise, ix_params_noise])$values
## [1] 3.389182e+01 1.856179e+01 3.248406e+00 2.474626e+00 1.958111e+00
## [6] 4.496991e-01 1.668400e-01 6.958517e-02 4.574457e-06 -4.585225e-07
eigen(hopt[ix_params_noise_b, ix_params_noise_b])$values
## [1] 33.891818 18.561787 3.248406 1.958111
eigen(hopt[ix_params_noise_sgt, ix_params_noise_sgt])$values
## [1] 2.474626e+00 4.496991e-01 1.668400e-01 6.958517e-02 4.574457e-06
## [6] -4.585225e-07
eigen(hopt[ix_params_noise_sgt[1:3], ix_params_noise_sgt[1:3]])$values
## [1] 2.474626e+00 1.668400e-01 -4.585225e-07
hopt[ix_params_noise_sgt[1:3], ix_params_noise_sgt[1:3]] %>% zapsmall()
## [,1] [,2] [,3]
## [1,] 2.4657285 -0.1430155 3e-07
## [2,] -0.1430155 0.1757371 2e-07
## [3,] 0.0000003 0.0000002 -5e-07
hopt[ix_params_noise_sgt[1:3], ix_params_noise_sgt[1:3]] %>% ginv() %>% zapsmall()
## [,1] [,2] [,3]
## [1,] 0.4 0.3 0.4
## [2,] 0.3 6.0 3.4
## [3,] 0.4 3.4 -2180917.9
eigen(hopt[ix_params_noise_sgt[4:6], ix_params_noise_sgt[4:6]])$values
## [1] 4.496991e-01 6.958517e-02 4.574457e-06
hopt[ix_params_noise_sgt[4:6], ix_params_noise_sgt[4:6]] %>% zapsmall()
## [,1] [,2] [,3]
## [1,] 0.4423478 -0.0523373 -0.0010726
## [2,] -0.0523373 0.0768619 0.0023817
## [3,] -0.0010726 0.0023817 0.0000791
hopt[ix_params_noise_sgt[4:6], ix_params_noise_sgt[4:6]] %>% ginv() %>% zapsmall()
## [,1] [,2] [,3]
## [1,] 2.86 11.08 -294.85
## [2,] 11.08 236.46 -6967.65
## [3,] -294.85 -6967.65 218382.47
par(mfcol = c(3,2))
opg[ix_params_noise, ix_params_noise] %>% corrplot(is.corr = FALSE)
ha[ix_params_noise, ix_params_noise] %>% corrplot(is.corr = FALSE)
hopt[ix_params_noise, ix_params_noise] %>% corrplot(is.corr = FALSE)
opg[ix_params_noise, ix_params_noise] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
ha[ix_params_noise, ix_params_noise] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
hopt[ix_params_noise, ix_params_noise] %>%cov2cor() %>% corrplot() %>% suppressWarnings()
B matrix
par(mfcol = c(3,2))
opg[ix_params_noise_b, ix_params_noise_b] %>% corrplot(is.corr = FALSE)
ha[ix_params_noise_b, ix_params_noise_b] %>% corrplot(is.corr = FALSE)
hopt[ix_params_noise_b, ix_params_noise_b] %>% corrplot(is.corr = FALSE)
opg[ix_params_noise_b, ix_params_noise_b] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
ha[ix_params_noise_b, ix_params_noise_b] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
hopt[ix_params_noise_b, ix_params_noise_b] %>%cov2cor() %>% corrplot() %>% suppressWarnings()
Distribution parameters
par(mfcol = c(3,2))
opg[ix_params_noise_sgt, ix_params_noise_sgt] %>% corrplot(is.corr = FALSE)
ha[ix_params_noise_sgt, ix_params_noise_sgt] %>% corrplot(is.corr = FALSE)
hopt[ix_params_noise_sgt, ix_params_noise_sgt] %>% corrplot(is.corr = FALSE)
opg[ix_params_noise_sgt, ix_params_noise_sgt] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
ha[ix_params_noise_sgt, ix_params_noise_sgt] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
hopt[ix_params_noise_sgt, ix_params_noise_sgt] %>%cov2cor() %>% corrplot() %>% suppressWarnings()
We compare the diagonal elements of
The columns starting with diff_ contain the absolute difference between the different objects
(
comp_diag = tibble(name = c(rep("AR", 4),
rep("MA_bwd", 4),
rep("MA_fwd", 4),
rep("B", 4),
rep("distr", 6)),
type = c(rep("system", 12),
rep("noise", 10)),
a_opg_diag = diag(opg),
b_ha_diag = diag(ha),
c_hopt_diag = diag(hopt))
)
comp_diag %>%
dtable()
## This version of bslib is designed to work with rmarkdown version 2.7 or higher.
We need the following helper function which calculates the sum of the absolute differences between the elements pertaining to (OPG, Hessian optim, Hessian analytic)
obtain_diff_abs_mat = function(df){
mat = df %>%
as.matrix()
names_df = colnames(df)
res = matrix(nrow = ncol(df), ncol = ncol(df))
for (ix_row in 1:ncol(df)){
for (ix_col in 1:ncol(df)) {
res[ix_row, ix_col] = (mat[, ix_row] - mat[, ix_col]) %>% abs() %>% sum(na.rm = TRUE)
}
}
rownames(res) = names_df
colnames(res) = names_df
return(res)
}
Investigate the differences between the diagonal elements for each group (AR, MA backwards, MA forwards, B, distribution).
comp_diag %>% select(contains("diag")) %>% obtain_diff_abs_mat()
## a_opg_diag b_ha_diag c_hopt_diag
## a_opg_diag 0.0000 46.801195 52.865005
## b_ha_diag 46.8012 0.000000 9.505489
## c_hopt_diag 52.8650 9.505489 0.000000
Plot AR parameters
comp_diag %>%
filter(name == "AR") %>%
rowid_to_column() %>%
select(-contains("diff")) %>%
select(-name, -type) %>%
pivot_longer(-rowid) %>%
ggplot(aes(x = rowid, y = value, fill = name)) +
geom_bar(stat = "identity", position = position_dodge()) +
ggtitle("Diagonal Elements for System Parameters: AR")
Plot MA backwards parameters
comp_diag %>%
filter(name == "MA_bwd") %>%
rowid_to_column() %>%
select(-contains("diff")) %>%
select(-name, -type) %>%
pivot_longer(-rowid) %>%
ggplot(aes(x = rowid, y = value, fill = name)) +
geom_bar(stat = "identity", position = position_dodge()) +
ggtitle("Diagonal Elements for System Parameters: MA backwards")
Plot MA forwards parameters
comp_diag %>%
filter(name == "MA_fwd") %>%
rowid_to_column() %>%
select(-contains("diff")) %>%
select(-name, -type) %>%
pivot_longer(-rowid) %>%
ggplot(aes(x = rowid, y = value, fill = name)) +
geom_bar(stat = "identity", position = position_dodge()) +
ggtitle("Diagonal Elements for System Parameters: MA forwards")
Plot noise parameters
comp_diag %>%
filter(name == "B") %>%
rowid_to_column() %>%
select(-contains("diff")) %>%
select(-name, -type) %>%
pivot_longer(-rowid) %>%
ggplot(aes(x = rowid, y = value, fill = name)) +
geom_bar(stat = "identity", position = position_dodge()) +
ggtitle("Diagonal Elements for Noise Parameters: B")
comp_diag %>%
filter(name == "distr") %>%
rowid_to_column() %>%
select(-contains("diff")) %>%
select(-name, -type) %>%
pivot_longer(-rowid) %>%
ggplot(aes(x = rowid, y = value, fill = name)) +
geom_bar(stat = "identity", position = position_dodge()) +
ggtitle("Diagonal Elements for Noise Parameters: Distribution")
hopt %>% diag() %>% sqrt()
## Warning in sqrt(.): NaNs produced
## [1] 2.608967522 4.923165175 2.623907143 8.288660201 2.619247677 4.845358745
## [7] 1.099944854 1.595651871 2.201818164 4.612020757 0.605369915 1.367756995
## [13] 1.629354911 4.393503760 1.936669761 5.652588739 1.570263821 0.419210130
## [19] NaN 0.665092310 0.277239802 0.008894932
Difference between all possible options for calculating the standard deviations of the parameters.
comp_sd = tibble(name = c(rep("AR", 4),
rep("MA_bwd", 4),
rep("MA_fwd", 4),
rep("B", 4),
rep("distr", 6)),
type = c(rep("system", 12),
rep("noise", 10))) %>%
unite("type_name", type, name, remove = FALSE) %>%
mutate(aa_opg_inv = diag(ginv(opg)) %>% sqrt(),
ab_hopt_inv = diag(ginv(hopt)) %>% sqrt(),
ac_ha_inv = diag(ginv(ha)) %>% sqrt(),
ba_rob_hopt = diag(ginv(hopt) %*% opg %*% ginv(hopt)) %>% sqrt(),
bb_rob_ha = diag(ginv(ha) %*% opg %*% ginv(ha)) %>% sqrt())
comp_sd %>%
mutate(across(where(is.numeric), ~ round(.x, digits = 3))) %>%
dtable()
comp_sd %>%
select(starts_with("a"), starts_with("b")) %>%
obtain_diff_abs_mat()
## aa_opg_inv ab_hopt_inv ac_ha_inv ba_rob_hopt bb_rob_ha
## aa_opg_inv 0.00000 70.489367 67.532362 249.89354 231.34195
## ab_hopt_inv 70.48937 0.000000 9.927695 182.59183 172.93476
## ac_ha_inv 67.53236 9.927695 0.000000 189.43229 166.10751
## ba_rob_hopt 249.89354 182.591831 189.432286 0.00000 31.76549
## bb_rob_ha 231.34195 172.934764 166.107507 31.76549 0.00000
comp_sd %>%
filter(type == "system") %>%
select(starts_with("a"), starts_with("b")) %>%
obtain_diff_abs_mat()
## aa_opg_inv ab_hopt_inv ac_ha_inv ba_rob_hopt bb_rob_ha
## aa_opg_inv 0.000000 8.175851 5.179494 32.95936 14.29749
## ab_hopt_inv 8.175851 0.000000 9.816698 27.93706 18.16038
## ac_ha_inv 5.179494 9.816698 0.000000 34.80916 11.37410
## ba_rob_hopt 32.959362 27.937064 34.809156 0.00000 31.57388
## bb_rob_ha 14.297493 18.160381 11.374097 31.57388 0.00000
comp_sd %>%
filter(type == "noise") %>%
select(starts_with("a"), starts_with("b")) %>%
obtain_diff_abs_mat()
## aa_opg_inv ab_hopt_inv ac_ha_inv ba_rob_hopt bb_rob_ha
## aa_opg_inv 0.00000 62.3135158 62.3528677 216.9341786 217.0444585
## ab_hopt_inv 62.31352 0.0000000 0.1109974 154.6547665 154.7743831
## ac_ha_inv 62.35287 0.1109974 0.0000000 154.6231304 154.7334103
## ba_rob_hopt 216.93418 154.6547665 154.6231304 0.0000000 0.1916113
## bb_rob_ha 217.04446 154.7743831 154.7334103 0.1916113 0.0000000
comp_sd %>%
filter(name == "B") %>%
select(starts_with("a"), starts_with("b")) %>%
obtain_diff_abs_mat()
## aa_opg_inv ab_hopt_inv ac_ha_inv ba_rob_hopt bb_rob_ha
## aa_opg_inv 0.0000000 0.59969563 0.60340252 0.85979233 0.85665695
## ab_hopt_inv 0.5996956 0.00000000 0.01339046 0.29419986 0.30039404
## ac_ha_inv 0.6034025 0.01339046 0.00000000 0.29820929 0.29507391
## ba_rob_hopt 0.8597923 0.29419986 0.29820929 0.00000000 0.02700318
## bb_rob_ha 0.8566569 0.30039404 0.29507391 0.02700318 0.00000000
comp_sd %>%
filter(name == "distr") %>%
select(starts_with("a"), starts_with("b")) %>%
obtain_diff_abs_mat()
## aa_opg_inv ab_hopt_inv ac_ha_inv ba_rob_hopt bb_rob_ha
## aa_opg_inv 0.00000 61.71382021 61.74946520 216.0743863 216.1878016
## ab_hopt_inv 61.71382 0.00000000 0.09760692 154.3605666 154.4739890
## ac_ha_inv 61.74947 0.09760692 0.00000000 154.3249211 154.4383364
## ba_rob_hopt 216.07439 154.36056662 154.32492111 0.0000000 0.1646081
## bb_rob_ha 216.18780 154.47398903 154.43833639 0.1646081 0.0000000
comp_sd %>%
filter(name == "AR") %>%
rowid_to_column() %>%
select(-contains("diff")) %>%
select(-name, -type, -type_name) %>%
pivot_longer(-rowid) %>%
ggplot(aes(x = rowid, y = value, fill = name)) +
geom_bar(stat = "identity", position = position_dodge()) +
ggtitle("Standard Errors of AR Parameters")
comp_sd %>%
filter(name == "MA_bwd") %>%
rowid_to_column() %>%
select(-contains("diff")) %>%
select(-name, -type, -type_name) %>%
pivot_longer(-rowid) %>%
ggplot(aes(x = rowid, y = value, fill = name)) +
geom_bar(stat = "identity", position = position_dodge()) +
ggtitle("Standard Errors of MA backwards Parameters")
comp_sd %>%
filter(name == "MA_fwd") %>%
rowid_to_column() %>%
select(-contains("diff")) %>%
select(-name, -type, -type_name) %>%
pivot_longer(-rowid) %>%
ggplot(aes(x = rowid, y = value, fill = name)) +
geom_bar(stat = "identity", position = position_dodge()) +
ggtitle("Standard Errors of MA forwards Parameters")
comp_sd %>%
filter(name == "B") %>%
rowid_to_column() %>%
select(-contains("diff")) %>%
select(-name, -type, -type_name) %>%
pivot_longer(-rowid) %>%
ggplot(aes(x = rowid, y = value, fill = name)) +
geom_bar(stat = "identity", position = position_dodge()) +
ggtitle("Standard Errors of B Parameters")
PARAMS[ix_params_noise_sgt]
## [1] -0.5243178 1.4924937 85616.3114036 0.1492405 1.9211332
## [6] 7.8907217
comp_sd %>%
filter(name == "distr") %>%
rowid_to_column() %>%
select(-contains("diff")) %>%
select(-name, -type, -type_name) %>%
pivot_longer(-rowid) %>%
ggplot(aes(x = rowid, y = value, fill = name)) +
geom_bar(stat = "identity", position = position_dodge()) +
ggtitle("Standard Errors of Distribution Parameters")
We choose the OPG for standard deviations in the correctly specified setting and the robust version calculated with analytic Hessian in the misspecified setting
opt_mod_params = fill_tmpl_whf_rev(theta = PARAMS,
tmpl = TMPL)
opt_mod_sd_correctlyspecified = fill_tmpl_whf_rev(theta = comp_sd %>% pull(aa_opg_inv),
tmpl = TMPL)
opt_mod_sd_misspecified = fill_tmpl_whf_rev(theta = comp_sd %>% pull(ba_rob_hopt),
tmpl = TMPL)
Print AR parameters with non-robust and robust standard errors in wide matrix format.
opt_mod_params$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 -0.107 -0.277
## [2,] 0 1 0.492 -1.013
opt_mod_sd_correctlyspecified$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 2.782 1.104
## [2,] 0 1 1.559 0.575
opt_mod_sd_misspecified$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 2.921 1.308
## [2,] 0 1 6.327 2.331
Print MA backwards parameters with non-robust and robust standard errors in wide matrix format.
opt_mod_params$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0.065 0.749
## [2,] 0 1 0.259 0.002
opt_mod_sd_correctlyspecified$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 2.615 2.985
## [2,] 0 1 1.751 1.057
opt_mod_sd_misspecified$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 3.689 10.528
## [2,] 0 1 7.864 3.656
Print MA forwards parameters with non-robust and robust standard errors in wide matrix format.
opt_mod_params$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0.561 1.671
## [2,] 0 1 -0.241 -0.296
opt_mod_sd_correctlyspecified$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0.985 2.617
## [2,] 0 1 0.328 0.932
opt_mod_sd_misspecified$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 3.203 7.564
## [2,] 0 1 0.780 2.077
Print static matrix B with non-robust and robust standard errors.
opt_mod_params$B %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] 0.936 0.289
## [2,] -0.223 0.182
opt_mod_sd_correctlyspecified$B %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] 0.693 1.088
## [2,] 0.279 0.260
opt_mod_sd_misspecified$B %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] 0.638 0.376
## [2,] 0.287 0.176
Print distribution parameters with non-robust and robust standard errors , \(\lambda, p, q\), for each variable (in columns)
opt_mod_params$distr %>% matrix(ncol = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] -0.524 0.149
## [2,] 1.492 1.921
## [3,] 85616.311 7.891
opt_mod_sd_correctlyspecified$distr %>% matrix(ncol = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] 1.094 1.419
## [2,] 2.826 12.576
## [3,] 0.000 409.497
opt_mod_sd_misspecified$distr %>% matrix(ncol = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] 0.444 2.295
## [2,] 2.261 21.203
## [3,] 0.000 614.852
Under correct specification may also choose, e.g., the analytic version of the Hessian to obtain the information matrix and calculate a robust version by using again the analytic version of the Hessian as well as the OPG
opt_mod_params = fill_tmpl_whf_rev(theta = PARAMS,
tmpl = TMPL)
opt_mod_sd_correctlyspecified_alt = fill_tmpl_whf_rev(theta = comp_sd %>% pull(ac_ha_inv),
tmpl = TMPL)
opt_mod_sd_misspecified_alt = fill_tmpl_whf_rev(theta = comp_sd %>% pull(bb_rob_ha),
tmpl = TMPL)
Print AR parameters with non-robust and robust standard errors in wide matrix format.
opt_mod_params$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 -0.107 -0.277
## [2,] 0 1 0.492 -1.013
opt_mod_sd_correctlyspecified_alt$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 2.626 1.294
## [2,] 0 1 0.847 0.289
opt_mod_sd_misspecified_alt$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 4.801 2.473
## [2,] 0 1 1.095 0.187
Print MA backwards parameters with non-robust and robust standard errors in wide matrix format.
opt_mod_params$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0.065 0.749
## [2,] 0 1 0.259 0.002
opt_mod_sd_correctlyspecified_alt$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 2.332 1.691
## [2,] 0 1 0.967 1.604
opt_mod_sd_misspecified_alt$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 3.933 1.605
## [2,] 0 1 1.409 3.492
Print MA forwards parameters with non-robust and robust standard errors in wide matrix format.
opt_mod_params$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0.561 1.671
## [2,] 0 1 -0.241 -0.296
opt_mod_sd_correctlyspecified_alt$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 1.055 3.111
## [2,] 0 1 0.353 1.272
opt_mod_sd_misspecified_alt$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 1.588 4.675
## [2,] 0 1 0.510 2.671
Print static matrix B with non-robust and robust standard errors.
opt_mod_params$B %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] 0.936 0.289
## [2,] -0.223 0.182
opt_mod_sd_correctlyspecified_alt$B %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] 0.644 0.616
## [2,] 0.258 0.198
opt_mod_sd_misspecified_alt$B %>% matrix(nrow = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] 0.644 0.364
## [2,] 0.279 0.177
Print distribution parameters with non-robust and robust standard errors , \(\lambda, p, q\), for each variable (in columns)
opt_mod_params$distr %>% matrix(ncol = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] -0.524 0.149
## [2,] 1.492 1.921
## [3,] 85616.311 7.891
opt_mod_sd_correctlyspecified_alt$distr %>% matrix(ncol = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] 0.683 1.690
## [2,] 2.379 15.377
## [3,] 0.000 467.316
opt_mod_sd_misspecified_alt$distr %>% matrix(ncol = 2) %>% round(digits = 3)
## [,1] [,2]
## [1,] 0.469 2.295
## [2,] 2.124 21.204
## [3,] 0.000 614.855
sgt_skewness = function(l, p, q){
v = q^(-1/p) *
(
(3*l^2+1) *
(beta(3/p,q-2/p)/beta(1/p,q)) -
4*l^2 *
(beta(2/p,q-1/p)/beta(1/p,q))^2
)^(-1/2)
skewness = 2 * q^(3/p) * l * v^3 / beta(1/p, q)^3 *
(8*l^2*beta(2/p,q-1/p)^3 -
3*(1+3*l^2) * beta(1/p,q) * beta(2/p, q-1/p) * beta(3/p,q-2/p) +
2 * (1+l^2) * beta(1/p,q)^2 * beta(4/p, q-3/p)
)
return(skewness)
}
sgt_kurtosis = function(l, p, q){
v = q^(-1/p) *
(
(3*l^2+1) *
(beta(3/p,q-2/p)/beta(1/p,q)) -
4*l^2 *
(beta(2/p,q-1/p)/beta(1/p,q))^2
)^(-1/2)
kurtosis = q^(4/p) * v^4 / beta(1/p, q)^4 *
(48*l^4*beta(2/p,q-1/p)^4 +
24 * l^2 * (1+3*l^2) * beta(1/p,q) * beta(2/p, q-1/p)^2 * beta(3/p,q-2/p) -
32 * l^2 * (1+l^2) * beta(1/p,q)^2 * beta(2/p, q-1/p) * beta(4/p,q-3/p) +
(1+10*l^2+5*l^4) * beta(1/p,q)^3 * beta(5/p, q-4/p)
)
return(kurtosis)
}
sgt_skewness(-0.52,1.49,85616)
## [1] -1.035358
sgt_kurtosis(-0.52,1.49,85616)
## [1] 6.318257
sgt_skewness(0.15,1.92,7.89)
## [1] 0.3251307
sgt_kurtosis(0.15,1.92,7.89)
## [1] 3.769704